Zadanie

Dopasuj model (lub modele) przewidujący wartość zmiennej $ConcreteStrength$ (Wytrzymałość na ściskanie betonu) . Model zbuduj w wybrany przez siebie sposób używając regresji wielorakiej. Oceń jakość modelu i spełnienie założeń. Porównaj otrzymane modele.

Rozwiązanie

Chociaż nie jesteśmy ekspertami w dziedzinie cementu i betonu, zamierzamy zastosować model regresji wielorakiej, by analizować i interpretować dane zawarte w załączonym pliku. Naszym celem jest zrozumienie, jak różne składniki wpływają na wytrzymałość betonu, wykorzystując dostępne dane.

Informacje o zbiorze danych

Zbiór danych “Concrete Compressive Strength” to zestaw informacji dotyczących składników i wytrzymałości betonu. Zawiera zmienne odnoszące się do składników mieszanki betonowej oraz ich wpływu na wytrzymałość betonu. Zawiera 9 zmiennych:

Cement (component 1)(kg in a m3 mixture): Ilość cementu w kilogramach na metr sześcienny mieszanki.
Blast Furnace Slag (component 2)(kg in a m3 mixture): Ilość żużla wielkopiecowego w kilogramach na metr sześcienny mieszanki.
Fly Ash (component 3)(kg in a m3 mixture): Ilość popiołu lotnego w kilogramach na metr sześcienny mieszanki.
Water (component 4)(kg in a m3 mixture): Ilość wody w kilogramach na metr sześcienny mieszanki.
Superplasticizer (component 5)(kg in a m3 mixture): Ilość superplastyfikatora w kilogramach na metr sześcienny mieszanki.
Coarse Aggregate (component 6)(kg in a m3 mixture): Ilość kruszywa grubego w kilogramach na metr sześcienny mieszanki.
Fine Aggregate (component 7)(kg in a m3 mixture): Ilość kruszywa drobnego w kilogramach na metr sześcienny mieszanki.
Age (day): Wiek mieszanki betonowej w dniach.
Concrete compressive strength(MPa, megapascals): Wytrzymałość na ściskanie betonu w megapaskalach.

Wczytanie zbioru danych

concrete_data <- read.csv("ConcreteData.csv", header=TRUE)
head(concrete_data)

Zmiana nazw zmiennych (na krótsze):

names(concrete_data) <- c("Cement", "BlastFurnaceSlag", "FlyAsh", "Water", "Superplasticizer", "CoarseAggregate", "FineAggregate", "Age", "ConcreteStrength")
head(concrete_data)

Następnie postaramy się dopasować model (lub modele), które bedą przewidywały w sposób liniowy wartość zmiennej $ConcreteStrength$. Aby to zrobić, potrzebujemy znaleźć zmienną objaśniające, które najbardziej odpowiadają do modelu regresji liniowej wielorakiej. Najpierw zapoznajmy się z naszymi danymi:

Ogólne statystyki:

summary(concrete_data)

##      Cement      BlastFurnaceSlag     FlyAsh           Water      
##  Min.   :102.0   Min.   :  0.00   Min.   :  0.00   Min.   :121.8  
##  1st Qu.:190.7   1st Qu.:  0.00   1st Qu.:  0.00   1st Qu.:166.6  
##  Median :255.0   Median : 19.00   Median :  0.00   Median :185.0  
##  Mean   :273.9   Mean   : 71.07   Mean   : 57.87   Mean   :181.8  
##  3rd Qu.:339.0   3rd Qu.:140.32   3rd Qu.:118.60   3rd Qu.:192.0  
##  Max.   :540.0   Max.   :316.10   Max.   :200.10   Max.   :247.0  
##  Superplasticizer CoarseAggregate  FineAggregate        Age        
##  Min.   : 0.000   Min.   : 801.0   Min.   :594.0   Min.   :  1.00  
##  1st Qu.: 0.000   1st Qu.: 932.0   1st Qu.:744.3   1st Qu.:  7.00  
##  Median : 6.700   Median : 968.0   Median :780.0   Median : 28.00  
##  Mean   : 6.141   Mean   : 974.4   Mean   :777.1   Mean   : 38.47  
##  3rd Qu.:10.100   3rd Qu.:1030.2   3rd Qu.:822.6   3rd Qu.: 28.00  
##  Max.   :32.200   Max.   :1145.0   Max.   :992.6   Max.   :365.00  
##  ConcreteStrength
##  Min.   : 6.27   
##  1st Qu.:23.20   
##  Median :33.42   
##  Mean   :34.98   
##  3rd Qu.:44.66   
##  Max.   :82.60

Dokładniejsze statystyki:

summary_data<-skim(concrete_data)
summary_data

Data summary
Name	concrete_data
Number of rows	908
Number of columns	9
_______________________
Column type frequency:
numeric	9
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Cement	1	273.93	99.23	102.00	190.70	255.00	339.00	540.0	▆▇▆▃▁
BlastFurnaceSlag	1	71.07	84.65	0.00	0.00	19.00	140.32	316.1	▇▂▂▁▁
FlyAsh	1	57.87	64.58	0.00	0.00	0.00	118.60	200.1	▇▁▃▂▁
Water	1	181.84	19.73	121.80	166.60	185.00	192.00	247.0	▁▅▇▂▁
Superplasticizer	1	6.14	5.50	0.00	0.00	6.70	10.10	32.2	▇▇▁▁▁
CoarseAggregate	1	974.41	76.56	801.00	932.00	968.00	1030.25	1145.0	▂▅▇▅▂
FineAggregate	1	777.12	74.02	594.00	744.27	780.05	822.65	992.6	▂▃▇▃▁
Age	1	38.47	47.54	1.00	7.00	28.00	28.00	365.0	▇▁▁▁▁
ConcreteStrength	1	34.98	16.14	6.27	23.20	33.42	44.66	82.6	▅▇▆▃▁

Histogramy dla poszczególnej zmiennej:

par(mfrow = c(3,4))

for (i in 1:ncol(concrete_data)){
  column_name<-colnames(concrete_data)[i]
  hist(concrete_data[[column_name]], 
       main = paste("Histogram - ", column_name), 
       xlab = "Value", ylab = "Frequecy")
}

Możemy przystąpić do wyznaczania zmiennych objaśniających do naszych modeli. Obliczymy współczynniki korelacji Pearsona $ConcreteStrength$ z każdą poszczególną zmienną w naszej bazie danych oraz stworzymy wykresy pokazujące każdą zmienną do zmiennej $ConcreteStrength$ (analiza wykresów i korelacji):

correlation_concrete <- cor(concrete_data$ConcreteStrength, concrete_data)
print(correlation_concrete)

##         Cement BlastFurnaceSlag      FlyAsh      Water Superplasticizer
## [1,] 0.4473283        0.1849669 -0.09264568 -0.2786437        0.3595313
##      CoarseAggregate FineAggregate       Age ConcreteStrength
## [1,]      -0.2056304    -0.1488172 0.3969586                1

for (col_name in names(concrete_data)) {
  p <- ggplot(data = concrete_data, aes_string(x = col_name, y = "ConcreteStrength")) +
    geom_point() + labs(x = col_name, y = "ConcreteStrength") + ggtitle(paste("ConcreteStrength względem", col_name," (Korelacja:", cor(concrete_data$ConcreteStrength, concrete_data[[col_name]]), ")"))
  print(p)
}

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Również stworzymy macierz korelacji, aby lepiej również zrozumieć jak wszystkie dane są ze sobą skorelowane:

ggcorrplot(cor(concrete_data), type='upper')

corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)

Najpierw stwórzmy pierwszy model, gdzie będą uwzględnione wszystkie zmienne:

model1 <- lm(ConcreteStrength ~ ., data = concrete_data)
summary(model1)

## 
## Call:
## lm(formula = ConcreteStrength ~ ., data = concrete_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.050  -5.071   0.829   5.674  26.235 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      10.627526  26.526229   0.401   0.6888    
## Cement            0.116777   0.008190  14.259  < 2e-16 ***
## BlastFurnaceSlag  0.097201   0.009827   9.892  < 2e-16 ***
## FlyAsh            0.075877   0.012028   6.309 4.41e-10 ***
## Water            -0.212400   0.041257  -5.148 3.23e-07 ***
## Superplasticizer  0.244278   0.096365   2.535   0.0114 *  
## CoarseAggregate   0.005622   0.009247   0.608   0.5434    
## FineAggregate     0.007508   0.010535   0.713   0.4762    
## Age               0.178619   0.006530  27.352  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.059 on 899 degrees of freedom
## Multiple R-squared:  0.6878, Adjusted R-squared:  0.685 
## F-statistic: 247.6 on 8 and 899 DF,  p-value: < 2.2e-16

Wnioski z pierwszego modelu

$R^2$ na poziomie $0.6878$, oznacza iż model wyjaśnia około 68,8% obserwacji wartości zmiennej $ConcreteStrength$.
Istnieją wartości p-value blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.
Największy wpływ na model mają $Cement$, $BlastFurnaceSlag$, $FlyAsh$, $Water$, $Superplasticizer$ oraz $Age$.
Wartość $p-value$ dla $Intercept$ większa niż $0.05$ oznacza, że przecięcie nie jest istotnie różne od zera. Co by onzaczało, że siła wytrzymałości betonu nie potrzebuje żadnych składników. Możliwe zmiany później.

Analiza założeń

Założenie 1. Liniowa zależność między zmienną objaśnianą, a objaśniającymi:

Możemy stwierdzić na podstawie $Adj. R^2=0.685$ oraz $p-value for F-statistic: < 2.2e-16$, który wskazuje, że model jest istotny statystycznie, że nie musimy wykluczyć zależności liniowej.

avPlots(model1)

Założenie 2. Zerowa średnia:

summary_model1_residuals <- skim(model1$residuals) %>% mutate(complete_rate= NULL) %>% mutate(n_missing= NULL) %>% mutate(numeric.hist = NULL) #niepotrzebne statystyki na ten moment
summary_model1_residuals

Data summary
Name	model1$residuals
Number of rows	908
Number of columns	1
_______________________
Column type frequency:
numeric	1
________________________
Group variables	None

Variable type: numeric

skim_variable	mean	sd	p0	p25	p50	p75	p100
data	0	9.02	-38.05	-5.07	0.83	5.67	26.24

Można powiedzieć, że średnia jest równa 0. Sprawdźmy to również za pomocą testu t-Studenta:

t.test(model1$residuals)

## 
##  One Sample t-test
## 
## data:  model1$residuals
## t = -1.3262e-15, df = 907, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.5874179  0.5874179
## sample estimates:
##     mean of x 
## -3.969378e-16

$P-value=1$ oznacza, że nie ma wystarczających dowodów, aby odrzucić hipotezę zerową, co sugeruje brak istotnych statystycznie różnic. Zatem wektor losowy reszt w modelu ma średnią zero.

Założenie 3. Liniowa niezależność zmiennych objaśniających:

corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)

Możemy zobaczyć z macierzy korelacji, że np. zmienna $Superplasticizer$ i $Water$ są bardzo mocno ze sobą skorelowane, zatem są od siebie zależne. Założenie to odrzuca nasz model i musimy postąpić z korektą.

#Tworzenie drugiego modelu

Spróbujemy ulepszyć model1, aby całościowo nasz nowy model był lepszy oraz spełniał założenia. Skoro wystąpił problem, że dane zmienne objaśniające w modelu są od siebie zależne, należy im się przyjrzeć i poczynić odpowiednie kroki:

model1 <- lm(ConcreteStrength ~ ., data = concrete_data)
summary(model1)

## 
## Call:
## lm(formula = ConcreteStrength ~ ., data = concrete_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.050  -5.071   0.829   5.674  26.235 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      10.627526  26.526229   0.401   0.6888    
## Cement            0.116777   0.008190  14.259  < 2e-16 ***
## BlastFurnaceSlag  0.097201   0.009827   9.892  < 2e-16 ***
## FlyAsh            0.075877   0.012028   6.309 4.41e-10 ***
## Water            -0.212400   0.041257  -5.148 3.23e-07 ***
## Superplasticizer  0.244278   0.096365   2.535   0.0114 *  
## CoarseAggregate   0.005622   0.009247   0.608   0.5434    
## FineAggregate     0.007508   0.010535   0.713   0.4762    
## Age               0.178619   0.006530  27.352  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.059 on 899 degrees of freedom
## Multiple R-squared:  0.6878, Adjusted R-squared:  0.685 
## F-statistic: 247.6 on 8 and 899 DF,  p-value: < 2.2e-16

Widzimy, że $Pr(>|t|)$ dla zmiennych $CoarseAggregate$ oraz $FineAggregate$ $>0.05$. Są potencjalnymi kandydatami do wyrzucenia z modelu.

corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)

Widzimy, że istnieje multikolinearności (dla zmiennej $Water$ i $Superplasticizer$ najbardziej). Zastosujemy użycie współczynnika VIF na naszym modelu. VIF jest współczynnikiem, który definiujemy dla poszczególnego współczynnika zmiennej objaśnianej. Powstaje jako iloraz wariancji współczynnika $\beta_i$ w pełnym modelu do jego wariancji w modelu, który ma tylko $x_i$ jako zmienną objaśniającą. Możemy go wyrazić wzorem \[ VIF(\beta_i) = \frac{1}{1 - R^2_{X_j|X_{-j}} }\] gdzie przez $R^2_{X_j|X_{-j}}$ oznaczamy współczynnik $R^2$ modelu regresji, gdzie zmienną objaśnianą jest $X_j$, a objaśniającymi są wszystkie pozostałe zmienne niezależne z modelu pierwotnego. Można zauważyć z postaci tego współczynnika, że VIF jest zawsze większy od 1, a coraz większe wartości VIF wskazują na coraz większy stopień kolinearności. Praktyka pokazuje, że zaalarmowani powinniśmy być, gdy VIF danej zmiennej przekracza 5, a poważnie martwić, gdy przekracza 10. Zobaczmy jak to się prezentuje u nas:

car::vif(model1)

##           Cement BlastFurnaceSlag           FlyAsh            Water 
##         7.298968         7.647387         6.668511         7.325227 
## Superplasticizer  CoarseAggregate    FineAggregate              Age 
##         3.106059         5.539066         6.719696         1.065311

Wnioski z wyniku:

Pomimo,że $Cement$, $BlastFurnaceSlag$ mają powyżej wartości 5, mają przeciętnie wysoką korelacje dodatnią ze zmienną $ConcreteStrength$.
Pomimo,że $Water$, $Flyash$ mają powyżej wartości 5, mają przeciętnie wysoką korelacje ujemną ze zmienną $ConcreteStrength$.
Dla zmiennych $CoarseAggregate$ oraz $FineAggregate$ również wyszły wartości powyżej 5. Z racji, że obie te zmienne najmniej pasowały do modelu, możemy spróbować stworzyć nowy model bez nich.

model2 <- lm(ConcreteStrength ~ . -CoarseAggregate-FineAggregate, data = concrete_data)
summary(model2)

## 
## Call:
## lm(formula = ConcreteStrength ~ . - CoarseAggregate - FineAggregate, 
##     data = concrete_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.483  -5.032   0.858   5.614  26.049 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      28.627889   4.136370   6.921 8.51e-12 ***
## Cement            0.111785   0.004109  27.202  < 2e-16 ***
## BlastFurnaceSlag  0.091174   0.004802  18.986  < 2e-16 ***
## FlyAsh            0.069215   0.007494   9.236  < 2e-16 ***
## Water            -0.236652   0.020964 -11.288  < 2e-16 ***
## Superplasticizer  0.228643   0.087622   2.609  0.00922 ** 
## Age               0.178599   0.006525  27.372  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.052 on 901 degrees of freedom
## Multiple R-squared:  0.6876, Adjusted R-squared:  0.6855 
## F-statistic: 330.6 on 6 and 901 DF,  p-value: < 2.2e-16

Wnioski z drugiego modelu

Widzimy, że pomimo tego, że model się nie polepszył ale również się nie pogorszył znacząco, $Adj. R^2$=0.6876 vs. 0.6878.
Możemy zaobserwować, że wszystkie zmienne mają silne powiązanie z model z wyjątkiem $Superplasticizer$. Wyeliminujemy ją oraz wtedy również multikorelacje z $Water$.
Pozostałe wartości $p-value$ są blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.
Statystyka F również potwierdza istnienie związku pomiędzy zmienną $ConcreteStrength$, a zmiennynmi objaśniającymi.

Stwórzmy kolejny model:

model3 <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+Age, data = concrete_data)
summary(model3)

## 
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh + 
##     Water + Age, data = concrete_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.524  -5.217   0.850   5.649  26.627 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      33.945116   3.611198    9.40   <2e-16 ***
## Cement            0.115895   0.003808   30.44   <2e-16 ***
## BlastFurnaceSlag  0.096952   0.004275   22.68   <2e-16 ***
## FlyAsh            0.080109   0.006243   12.83   <2e-16 ***
## Water            -0.270327   0.016575  -16.31   <2e-16 ***
## Age               0.179726   0.006532   27.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.081 on 902 degrees of freedom
## Multiple R-squared:  0.6853, Adjusted R-squared:  0.6835 
## F-statistic: 392.8 on 5 and 902 DF,  p-value: < 2.2e-16

Również możemy zauważyć,że $R^2$ nie zmienił się zbytnio w stosunku do poprzedniego modelu.

Wnioski z modelu trzeciego

$R^2$ na poziomie $0.6853$, oznacza iż model wyjaśnia około 68,5% obserwacji wartości zmiennej $ConcreteStrength$.
Wartości p-value są blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.
Statystyka F również potwierdza istnienie związku pomiędzy zmienną $ConcreteStrength$, a zmiennynmi objaśniającymi.

Zobaczmy współczynnik VIF:

car::vif(model3)

##           Cement BlastFurnaceSlag           FlyAsh            Water 
##         1.570321         1.440448         1.788143         1.176643 
##              Age 
##         1.060627

Współczynnik ten jest poniżej 5 dla każdej zmiennej, zatem zmienne objaśniające nie są ze sobą skorelowane.

Zobaczmy, czy możemy ulepszyć powyższy model. Spójrzmy na histogramy naszych zmiennych:

par(mfrow = c(2,3))
names <- c("Cement","BlastFurnaceSlag","FlyAsh","Water","Age","ConcreteStrength")
for (i in names){
  hist(concrete_data[[i]], 
       main = paste("Histogram - ", i), 
       xlab = "Value", ylab = "Frequecy")
}

Możemy potraktować zmienną $Age$ logarytmem ze względu na jej wykładniczy charakter (zmienne $FlyAsh$ oraz $BlastFurnaceSlag$ mają wartości = 0).

model4 <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+log(Age), data = concrete_data)
summary(model4)

## 
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh + 
##     Water + log(Age), data = concrete_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.6524  -4.1841  -0.1371   4.2187  21.4391 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      10.551137   2.443128   4.319 1.74e-05 ***
## Cement            0.119674   0.002559  46.772  < 2e-16 ***
## BlastFurnaceSlag  0.090130   0.002865  31.456  < 2e-16 ***
## FlyAsh            0.063087   0.004183  15.081  < 2e-16 ***
## Water            -0.264682   0.011029 -23.999  < 2e-16 ***
## log(Age)          9.633018   0.182725  52.719  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.096 on 902 degrees of freedom
## Multiple R-squared:  0.8581, Adjusted R-squared:  0.8574 
## F-statistic:  1091 on 5 and 902 DF,  p-value: < 2.2e-16

Wnioski z czwartego modelu

$R^2$ (oraz $Adj. R^2$) znacząco podwyższyły się w porównaniu do poprzedniego modelu ($R^2$=0,6853 vs 0.8581)
$R^2$ na poziomie $0.8581$, oznacza iż model wyjaśnia około 85,8% obserwacji wartości zmiennej $ConcreteStrength$.
Wartości p-value są blisikie zeru, co pozwala nam stwiedzić iż istnieje związek pomiędzy zmiennymi.
Statystyka F również potwierdza istnienie związku pomiędzy zmienną $ConcreteStrength$, a zmiennynmi objaśniającymi.

Przeprowadźmy jeszcze raz analizę założeń modelu:

Druga analiza założeń

Założenie 1. Liniowa zależność między zmienną objaśnianą, a objaśniającymi:

Możemy stwierdzić na podstawie $Adj. R^2=0.8574$ oraz $p-value for F-statistic: < 2.2e-16$, który wskazuje, że model jest istotny statystycznie, że nie musimy wykluczyć zależności liniowej.

avPlots(model4)

Założenie 2. Zerowa średnia:

summary_model4_residuals <- skim(model4$residuals) %>% mutate(complete_rate= NULL) %>% mutate(n_missing= NULL) %>% mutate(numeric.hist = NULL)
summary_model4_residuals

Data summary
Name	model4$residuals
Number of rows	908
Number of columns	1
_______________________
Column type frequency:
numeric	1
________________________
Group variables	None

Variable type: numeric

skim_variable	mean	sd	p0	p25	p50	p75	p100
data	0	6.08	-16.65	-4.18	-0.14	4.22	21.44

Można powiedzieć, że średnia jest równa 0. Sprawdźmy to również za pomocą testu t-Studenta:

t.test(model4$residuals)

## 
##  One Sample t-test
## 
## data:  model4$residuals
## t = 2.9355e-15, df = 907, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3959577  0.3959577
## sample estimates:
##    mean of x 
## 5.922423e-16

$P-value=1$ oznacza, że nie ma wystarczających dowodów, aby odrzucić hipotezę zerową, co sugeruje brak istotnych statystycznie różnic. Zatem wektor losowy reszt w modelu ma średnią zero.

Założenie 3. Liniowo niezależne zmienne objaśniające:

corrplot(cor(concrete_data), method = "number", type = "upper", number.cex = 0.8, tl.cex = 0.8, tl.col = "black", diag = FALSE)

car::vif(model4)

##           Cement BlastFurnaceSlag           FlyAsh            Water 
##         1.573182         1.435734         1.781191         1.155967 
##         log(Age) 
##         1.038097

Współczynnik ten jest poniżej 5 dla każdej zmiennej, zatem zmienne objaśniające nie są ze sobą skorelowane.

Założenie 4. Homoskedastyczność:

plot(model4)

bptest(model4)

## 
##  studentized Breusch-Pagan test
## 
## data:  model4
## BP = 10.529, df = 5, p-value = 0.06156

Otrzymaliśmy $p-value$ znacznie większe niż 0.05, zatem nie mamy wystarczająco dowodów aby odrzucić hipotezę o homoskedastycznośći reszt.

Założenie 5. Normalność reszt:

Zbadamy to za pomocą histogramu reszt, wykresy kwantyl-kwantyl oraz testu Shapiro-Wilk’a:

ggplot(model4, aes(x=.resid)) + geom_histogram(bins=40) +labs(title = "Histogram reszt w modelu", x = "Wartość", y = "Częstotliwość")

ggplot(model4, aes(sample =.resid)) + geom_qq() + geom_qq_line() + labs(title = "Wykres kwantyl-kwantyl reszt", x = "Kwantyle teoretyczne", y = "Kwantyle próbkowe")

Histogram, oraz wykres kwantyl-kwantyl sugerują iż rozkład reszt może być normalny. Zweryfikujemy to za pomocą testu Shapiro-Wilk’a:

shapiro.test(model4$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model4$residuals
## W = 0.99769, p-value = 0.237

Otrzymaliśmy $p-value>0.05$, zatem nie mamy podstaw aby odrzucić hipotezę o normalności reszt.

Podsumowując, reszty w modelu spełniają klasyczne założenia modelu regresji.

Porównanie modeli

Porównajmy nasz model ostateczny z poprzednimi:

st_err1 = summary(model1)$sigma
st_err2 = summary(model2)$sigma
st_err3 = summary(model3)$sigma
st_err4 = summary(model4)$sigma
paste("Residual Standard Error :", "Model 1:",round(st_err1,4),"|", "Model 2:",round(st_err2,4),"|", "Model 3:",round(st_err3,4), "Model 4:",round(st_err4,4) )

## [1] "Residual Standard Error : Model 1: 9.0591 | Model 2: 9.0516 | Model 3: 9.0807 Model 4: 6.0963"

Model o mniejszym błędzie standardowym jest lepszy od innych, czyli u nas model ostateczny - model4.

Dodatkowo jednym z czynników porównującym model może być również współczynnik determinacji $Adj. R^2$, który tłumaczy jaki model lepiej wyjaśnia zmienność danych:

rsqr1 = summary(model1)$adj.r.squared
rsqr2 = summary(model2)$adj.r.squared
rsqr3 = summary(model3)$adj.r.squared
rsqr4 = summary(model4)$adj.r.squared
paste("Adjusted R^2 :", "Model 1:",round(rsqr1,4),"|", "Model 2:",round(rsqr2,4),"|", "Model 3:",round(rsqr3,4), "Model 4:",round(rsqr4,4) )

## [1] "Adjusted R^2 : Model 1: 0.685 | Model 2: 0.6855 | Model 3: 0.6835 Model 4: 0.8574"

Porównując współczynniki determinacji również możemy stwierdzić, że model4 najlepiej wyjaśnia zmienność danych.

Trenowanie modelu

Spróbujemy dodatko podzielić zbiór naszych danych na zbiór treningowy oraz zbiór testowy. Następnie porównamy nasz model z tym:

set.seed(1311)
split_data <- createDataPartition(concrete_data$ConcreteStrength, p = 0.75, list = FALSE)
train_data <- concrete_data[split_data, ]
test_data <- concrete_data[-split_data, ]

model4train <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+log(Age), data = train_data)
summary(model4train)

## 
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh + 
##     Water + log(Age), data = train_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.6395  -4.2858  -0.0804   4.3773  20.9613 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       8.537399   2.790595   3.059  0.00231 ** 
## Cement            0.119369   0.002933  40.695  < 2e-16 ***
## BlastFurnaceSlag  0.089144   0.003391  26.285  < 2e-16 ***
## FlyAsh            0.063314   0.004836  13.093  < 2e-16 ***
## Water            -0.251251   0.012638 -19.880  < 2e-16 ***
## log(Age)          9.519809   0.210946  45.129  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.138 on 677 degrees of freedom
## Multiple R-squared:  0.8553, Adjusted R-squared:  0.8543 
## F-statistic: 800.5 on 5 and 677 DF,  p-value: < 2.2e-16

Dodatkowym sposobem na poprawienie jakości modelu jest jeszcze spojrzeć na odległośći Cooka oraz dźwignie:

ggplot(model4train, aes(.hat, .stdresid))+ geom_point(aes(size=.cooksd))+ stat_smooth(method='loess', formula=y~x, se=FALSE)+ labs(title='Leverage vs Standardized Residuals',x='Leverage', y='Standardized Residuals', size='Cooks distance')

Postaramy się usunąć obserwacje o wysokiej dźwigni, oraz odległości Cooka:

train_data2 <- train_data%>% mutate(cooks_D = cooks.distance(model4train))%>% filter(cooks_D <= 3*mean(cooks_D))

Sprawdzimy czy nasz model na zmienionych danych:

model4train2 <- lm(ConcreteStrength ~ Cement+BlastFurnaceSlag+FlyAsh+Water+log(Age), data = train_data2)
summary(model4train2)

## 
## Call:
## lm(formula = ConcreteStrength ~ Cement + BlastFurnaceSlag + FlyAsh + 
##     Water + log(Age), data = train_data2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.8816  -3.8364  -0.1662   3.7222  15.4607 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      10.650197   2.580683   4.127 4.18e-05 ***
## Cement            0.118353   0.002617  45.231  < 2e-16 ***
## BlastFurnaceSlag  0.087462   0.003046  28.714  < 2e-16 ***
## FlyAsh            0.068838   0.004347  15.834  < 2e-16 ***
## Water            -0.259741   0.011980 -21.680  < 2e-16 ***
## log(Age)          9.373795   0.187549  49.980  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.196 on 619 degrees of freedom
## Multiple R-squared:  0.8899, Adjusted R-squared:  0.889 
## F-statistic:  1000 on 5 and 619 DF,  p-value: < 2.2e-16

#Ostateczne porównanie modeli:

st_err4 = summary(model4)$sigma
st_err4train = summary(model4train2)$sigma
paste("Residual Standard Error :", "Model 4:",round(st_err4,4),"|", "Model 4 treningowy:",round(st_err4train,4))

## [1] "Residual Standard Error : Model 4: 6.0963 | Model 4 treningowy: 5.1959"

Model o mniejszym błędzie standardowym jest lepszy, czyli u nas model4train2.

rsqr4 = summary(model4)$adj.r.squared
rsqr4train = summary(model4train2)$adj.r.squared
paste("Adjusted R-squared :", "Model 4:",round(rsqr4,4),"|", "Model 4 treningowy:",round(rsqr4train,4))

## [1] "Adjusted R-squared : Model 4: 0.8574 | Model 4 treningowy: 0.889"

Również współczynnik determinancji jest wyższy w modelu model4train2, zatem uwzględnienie odległości Cooka polepszyło nasz model, oraz jego parametry na zbiorze testowym.

Podsumowanie

Najlepsze dopasowanie naszego modelu, przy jednoczesnym spełnieniu założeń regresji wielorakiej otrzymaliśmy w modelu ze zmiennymi objaśniającymi: $Cement$, $BlastFurnaceSlag$, $FlyAsh$, $Water$ oraz $log(Age)$.

Dla tego samego modelu, poprzez nauczenie na zbiorze treningowym oraz uwzględnienie odległości Cooka otrzymaliśmy najlepsze dopasowanie, oraz najniższy bład standardowy na zbiorze testowym. Nasz model końcowy wyjaśnia około 88,9% zmiennej $ConcreteStrength$ w zbiorze testowym, co możemy uznać za bardzo dobre dopasowanie.

Projekt Regresja i Analiza Wariancji

Adrian Pilarczyk i Maria Bul

Zadanie

Rozwiązanie

Informacje o zbiorze danych

Wczytanie zbioru danych

Wnioski z pierwszego modelu

Analiza założeń

Wnioski z drugiego modelu

Wnioski z modelu trzeciego

Wnioski z czwartego modelu

Druga analiza założeń

Porównanie modeli

Trenowanie modelu

Podsumowanie